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ABSTRACT 

We study the field millisecond pulsar population to infer its intrinsic distribution in 
spin period and luminosity and to determine its spatial distribution within the Galaxy. 
Our likelihood analysis on data from extant surveys (22 pulsars with periods < 20 ms) 
accounts for the following important selection effects: (1) the survey sensitivity as a 
function of direction, spin period, and sky coverage; (2) interstellar scintillation, which 
modulates the pulsed flux and causes a net increase in search volume ~ 30%; and (3) 
errors in the pulsar distance scale. 

Adopting power-law models (with cutoffs) for the intrinsic distributions, the 
analysis yields a minimum period cutoff Pmin > 0.65 ms (99% conhdence), a period 
distribution oc p-2-0±o.33 pseudo-luminosity distribution oc (where 

Lp = flux density x distance^, for Lp > 1.1 mJy kpc^). 

We hnd that the column density of millisecond pulsars (uncorrected for beaming 
effects) is ~ 50^20 kpc“^ in the vicinity of the solar system. For a Gaussian model the 
z scale height is O.GS^o'}^ corresponding to local number density 291}]” kpc“^. 
(For an exponential model the scale height becomes kpc and the number 

density 44^^g kpc“^.) Estimates of the total number of MSPs in the disk of the Galaxy 
and for the associated birthrate are given. The contribution of a diffuse halo-like 
component (tracing the Galactic spheroid, the halo or the globular cluster density 
profile) to the local number density of MSPs is limited to ;^1% of the midplane value. 

We consider a kinematic model for the MSP spatial distribution in which objects 
in the disk are kicked once at birth and then orbit in a smooth Galactic potential, 
becoming dynamically well-mixed. The analysis yields a column density 49^^^ kpC-2 
(comparable to the above), a birth 2 kick velocity 52^|[ km s“^ and a 3D velocity 
dispersion of ~ 84 km s“^. MSP velocities are smaller than those of young, long-period 
pulsars by about a factor of 5. The kinematic properties of the MSP population 
are discussed, including expected transverse motions, the occurrence of asymmetric 
drift, the shape of the velocity ellipsoid and the z scale height at birth. If MSPs are 
long-lived then a significant contribution to observed MSP z velocities owes to diffusive 
processes that increase the scale height of old stellar populations; our best estimate of 
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the ID velocity kick that is unique to MSP evolution is ~ 40 km s“^ if such diffusion 
is taken into account. 

The scale heights of millisecond pulsars and low-mass X-ray binaries are consistent, 
suggesting a common origin and that the primary channel for forming both classes 
of objects imparts only low velocities. Binaries involving a common envelope phase 
and a neutron-star forming supernova explosion can yield such objects, even with 
explosion asymmetries like those needed to provide the velocity distribution of isolated, 
nonspunup radio pulsars. 

Future searches for MSPs may be optimized using the model results. As an 
example, we give the expected number of detectable MSPs per beam area and the 
volumes of the Galaxy sampled per beam area for a hypothetical Green Bank Telescope 
all sky survey. Estimates for the volume that must be surveyed to find a pulsar faster 
than 1.5 ms are given. We also briefly discuss how selection effects associated with fast 
binaries influence our results. 

Subject headings: pulsars, stars-binary: 


INTRODUCTION 


Millisecond pulsars (MSPs) differ from slower-spin pulsars in important ways. First, their 
spindown rates and derived surface magnetic fields are several orders of magnitude smaller. MSPs 
have implied fields of 10^'® — 10® Gauss, while pulsars with periods of order one second are 
characterized by magnetic fields of 10^^ — 10^^ Gauss. Closely related is the observation that the 
characteristic spin-down times of MSPs, ranging from several tenths to tens of Gyr, far exceed 
those of slower-spin pulsars. Some MSPs were born with periods near their present-day values 
and are, consequently, much younger than their spindown times (pamilo, Thorsett Sz Kulkarni 


1994). However, MSPs are thought to be active as radio pulsars for hundreds to thousands of 
times longer than strong-field pulsars. Taking active lifetimes into account, it appears that there 
may be comparable numbers of young pulsars and MSPs in the Galaxy though the birth rate of 
MSPs is ~ 10^ times smaller. 

A second significant difference is that more than 2/3 of MSPs are in binary systems, while 
young, strong-field pulsars are largely solitary objects, a fact which has both theoretical and 
observational implications. Clearly, the evolutionary pathways that gives rise to MSPs are 


integrally related to the interaction of the binary stars (Alpar et al. 1982, Ruderman &: Shaham 


1983 ; for a general review see Bhattacharya &: van den Heuvel 1991 ). Moreover, searches for MSPs 
must confront the additional selection effects that mitigate against the detection of accelerated 
pulsars. As is shown below, detection of the fastest spinning pulsars is inhibited by any effects that 
smear out the pulse, such as dispersive propagation in the interstellar medium. Orbital motion, 
uncompensated for in the surveys we analyze here, also smears out pulses according to the change 
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in velocity over the duration of the observation and is therefore most important for short period 
pulsars in compact binaries (e.g. pohnston &: Kulkarni 19^ ). 

In this paper, we analyze the spatial distribution of MSPs. Our purpose is to derive the best 
estimates for MSP population parameters through careful consideration of survey sensitivities as 
a function of pulse period, dispersion measure and other relevant factors. Our census essentially 
measures the local number density of MSPs, the fall off in number density above the Galactic 
plane and establishes upper limits on a diffuse, halo-like component of the MSP population. We 
analyze the implications of the spatial distribution for the kinematics of the MSP population, 
inferring the diffusive and impulsive velocity increments suffered. We compare predictions of the 
distribution of proper motions to extant observations, compare the spatial distribution of LMXBs 
and MSPs, and describe the importance of our determination of the MSP kick velocity for binary 
evolutionary scenarios. We provide detailed analysis of the influence of selection effects on the 
discovery of short period pulsars and pulsars in binaries. 

Our results have immediate relevance in a number of respects. Recent MSP surveys have 
been conducted on the premise that the MSPs are essentially isotropically distributed around the 
Sun, at least to the depths that surveys probe. Our results establish the scale height and show 
that Arecibo type surveys see beyond it. 


The third way that MSPs differ from slow-spin pulsars is in their peculiar space motions. 
This is one of the main conclusions of the present paper. Kinematic evidence (e.g. Dewey, Cordes 
& Wolszczan 1988; Cordes et al. 1990 jWolszczan 1994 ; Nice & Taylor 1995; Nicastro k, Johnston 
1995) suggests that MSPs are low velocity objects, with typical transverse speeds <100 — 200 km 
s“^ . Such velocities are much less than young pulsars, which have an average speed ~ 500 km 
s“^ (Lyne k Lorimer 1994; Cordes k Chernoff 1996). Our results allow us to put more stringent, 
albeit statistical, limits on the MSP velocities than has hitherto been achieved. 


The use of the spatial distribution of MSPs as an indirect means for determining their peculiar 
velocities is more robust than an analysis of proper motion data. The primary reason is that the 
orbits of MSPs are perturbed significantly from circular motion around the galactic center so that, 
given their ages (> 0.1 Gyr), corrections for differential galactic rotation cannot be made. Arnaud 
k Rothenflug (1981) applied a similar spatial analysis to young, high-held pulsars as a means for 
determining their velocities. (The methodology is correct but their assumption that high-held 
pulsars form a steady, relaxed population is not. Today we know that ~ 25-30% escape the Galaxy 
and that many radio pulsars shut off before traveling to the limiting distance for detection.) 


Our approach differs in several ways from those taken by other authors. First, we use a 
likelihood analysis to provide the best estimates of the MSP population parameters, to account 
accurately for survey selection effects and distance errors, and to express clearly our physical 
assumptions. Second, we restrict our analysis to pulsars with spin periods < 20 ms. We do 
so because it appears that these neutron stars (NS) are distinct from other pulsars that may 
have undergone accretion-driven spinup but were left with longer periods (Bhattacharya k van 
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den Heuvel 1991). We also consider them to be distinct from the higher mass NS-NS binaries, 
which have pulsars with longer pulse periods (B1913+16, P = 59 ms [Taylor & Weisberg 1989]; 
B1534+12, P = 38 ms [Wolszczan 1991]). Third, our approach includes the effects of interstellar 
scintillations, which modulate the pulsar flux density and influence the rate of detection in surveys 


The paper is organized as follows. In §|^ we discuss sensitivities of pulsar surveys and derive 
quantities that are needed in our analysis. A preliminary attack on the problem is given in 
where we present a ^VfVmax' analysis, analogous to that used on quasars and gamma-ray 
burst sources, in order to illustrate the uncertainties involved in pulsar surveys. We derive the 
survey likelihood function in and apply it in §P to a disk-only distribution of MSPs. In we 
apply the analysis to a disk model based on numerical integration of NS orbits in the galactic 
potential. We consider a combined disk and diffuse halo-like model in Sections |8[jII| present 
the implications of our results for the MSP birthrate, the origins of MSPs and their velocities, and 
for the optimization of MSP surveys. In §14 we summarize the paper. 


2. PULSAR SURVEYS 


2.1. Minimum Detectable Flux Smin 


A pulsar survey has a minimum detectable flux density Smin that depends on radiometer 
noise, the pulse shape, and details of the Fourier analysis used to find pulsars. In Appendix ^ we 
derive Smin for surveys, which includes pulse broadening effects (from interstellar dispersion and 
scattering and from detector time constants) and the effects of flux variations from interstellar 
scintillations (Appendix P), which are strong for the MSP surveys. 

For a given direction, the minimum detectable flux density is a function of pulse width and 
period and radiometer-noise level. The minimum detectable flux density is a function of direction 
(galactic coordinates i,b), dispersion measure [DM = dxne{x), where rig is the free-electron 
density and D is the pulsar distance], radio frequency {u), bandwidth (Ai^) and the number of 
channels (Nch), as well as system temperature (T^^^), telescope gain (G), the intrinsic pulse duty 
cycle, and additional survey dependent factors: 


Smin — Smin {£, b, P, DM, V, Au, Nch, Tsys, G,...). 


( 1 ) 


Smin depends on additional, unspecified parameters, especially those that describe orbital motion. 
In most of this paper we ignore such motion; however in we discuss survey biases against short 
period binaries and their possible effects on our conclusions. In Figure || we show Smin plotted 
against period for several values of DM. These curves apply for drift-scan surveys made with 
the Arecibo telescope at 0.43 GHz toward the galactic pole and at zero zenith angle. We have 
used a numerical version of a 408 MHz survey ( Haslam et al. 1982| ) to calculate the background 
sky temperature. We assume a Gaussian pulse shape with 3% intrinsic duty cycle. While this 
duty cycle is shorter than those of some MSPs, the duty cycle and, hence, Smin-, is dominated 
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by extrinsic pulse broadening from dispersion, scattering and instrumentation. Pulse broadening 
from scattering is taken into account by estimating it from large scale galactic models for the 
electron density, as described by Cordes et al. (1991) and Taylor &: Cordes (1993; hereafter TC). 

The maximum distance to which a particular pulsar is detectable, Dmax, is given by 


D 


max 


Lr 


-|l/2 


Sr, 


( 2 ) 


We use the ‘pseudo-luminosity’, Lp = SD^, that is often adopted in population studies of pulsars. 
Though it is preferable to use a physical luminosity (i.e. expressed in units of erg s“^) in analyzing 
pulsar statistics (Chernoff &: Cordes 1996a), estimation of physical luminosities for MSPs is not 
yet possible because we do not understand radio beaming in MSPs to the same extent that we 
do for young, strong-held objects (Backer 1976; Rankin 1983; Lyne & Manchester 1988; Rankin 
1993). Note that Smin depends implicitly on the dispersion measure in a given direction which, in 
turn, depends on Dmax- For this reason, Dmax must be found iteratively. 

Figure shows (as solid lines) Smin plotted against distance for several values of pulse period 
and for several directions. These curves illustrate the strong dependence of Smin on direction and 
period. The dashed line indicates the inverse square-law variation of flux density for a source of 
30 mJy at a distance of 1 kpc. The intersection points of the dashed and solid curves determine 
the maximum distances Dmax that a pulsar can be seen for the different cases. 


From Dmax, the total volume to which the survey is sensitive in a given beam area is 

1 3 

^max = '^^bDmax- 


( 3 ) 


Because of the strong period and luminosity dependence of Dmax and, hence, Vmax, the latter 
quantity may be used to determine the period and luminosity distributions of MSPs. 


2.2. MSP Surveys: Properties, Volumes and Distances Sampled 


We have applied the likelihood analysis to 8 pulsar surveys that have been reported in 
the literature. Six use the Arecibo telescope, yielding 11 MSPs, the seventh is the Parkes 
southern-hemisphere survey that yielded 10 MSPs, as reported by Manchester (1994) when the 
survey was about 75% complete. [| The eighth survey is the Jodrell Bank survey, a portion of 
which has been reported in the discovery of the MSP J1012-I-5307 (Nicastro et al. 1995). We 
have not included J0218-I-4232 ( Navarro et al. 1995| ) because it was discovered in an aperture 
synthesis survey with selection criteria quite different from the periodicity searches of the surveys 


^The final tally of the Parkes survey is 17 MSPs with all data analyzed, of which ~ 95% was relatively free of RFI (M. Bailes, private 
communication). 
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PULSE PERIOD (ms) 

Fig. 1.— Plot of minimum detectable flux density vs. pulse period for a typical survey at Arecibo. 
The different curves are for dispersion measures DM = 0 (dashed line) and DM = 10 x 2"', with 
n = 0,1,... (solid lines, left to right) and apply to observations toward the galactic poles and at 
zero zenith angle. 
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Fig. 2.— (Solid lines:) Minimum detectable flux density for an Arecibo survey for the labelled 
pulse periods (ms) and galactic latitudes (deg). The plotted values apply for sources viewed at the 
telescope beam center. At the beam half-power point, the values are doubled. (Dashed line;) The 
inverse-square law variation of flux density for a source with 30 mJy when at a 1-kpc distance. The 
distances where the solid lines cross the dashed line are values of D^ax, the maximum distance 
at which the particular object could be detected. At the beam half-power point, the maximum 
distances are smaller. 
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we consider. In a future analysis of the distributions of MSP orbital parameters, we will extend 
our analysis to synthesis surveys. 

Table ^ lists the surveys we have used. The columns include the (1) survey number (an 
arbitrary choice based solely on the order in which we analyzed the surveys); (2) observatory site 
for the survey; (3) survey frequency; (4) solid angle of interference-free observations; (5) number 
of MSPs found in the survey; (6) system temperature of the telescope, expressed in Janskys, for 
observations at the zenith and toward the Galactic poles; (7) minimum flux density for detection of 
long-period pulsars at the zenith and toward the Galactic poles; and (8) survey reference number. 

Table |2| lists the MSPs used in our analysis. Figure ^ shows the total volume searched in 
each of the 8 programs as a function of spin period for a fixed luminosity of 16 mJy kpc^. The 
Parkes survey by itself has searched the largest volume. The Arecibo surveys in aggregate cover a 
comparable volume, a few kpc^ for a period of 5 ms. For comparison, the lower panel in Figure 
shows the maximum distance surveyed, (Dmax), averaged over all directions searched in a given 
survey. The Arecibo surveys probe about 3 times more deeply than the Parkes survey. 

As we demonstrate in this paper, deep, high-latitude surveys at Arecibo (see references 
below) sample distances that are well beyond the scale height of the population, while the 
southern-hemisphere Parkes survey, covering a larger area on the sky but being less sensitive at 
most periods, is better optimized to finding MSPs. We note that the Arecibo search volumes that 
are devoid of MSPs provide the most stringent constraints on the scale height of MSPs, whereas 
the totality of detected MSPs essentially determines the local MSP number density. 


3. v/Vmax for millisecond pulsars 

The main method of analysis in this paper uses a likelihood function to determine intrinsic 
properties of the MSP population after accounting for survey selection effects embodied in Smin, 
as calculated above. Here, we motivate our discussion by applying a VjVmax analysis to MSP 
surveys. Consider the line of sight to a MSP discovered in a survey. Subsequent observations yield 
precise determinations of P, DM, £, b. The flux density is Sd at the time of discovery and S as 
a long time average. The flux density is time dependent, owing predominantly to refractive and 
diffractive interstellar scintillation (DISS; e.g. Rickett 1990; Kaspi &: Stinebring 1992, Stinebring 
&: Condon 1990; Cordes, Weisberg & Boriakoff 1985). A distance estimate derives from DM 
and the TC model for the interstellar electron density. The model and, in some cases, auxiliary 
measurements (timing parallax, neutral-hydrogen absorption, and association with supernova 
remnants; Frail & Weisberg 1990) yield a range of possible distances, \Dl,Du]. 

The volume between us and a given pulsar in a beam of solid angle is 

U = \^bD\ 


(4) 
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1 10 10 ® 
PERIOD (ms) 


Fig. 3.— (Top:) Volume searched as a function of spin period for a pseudo luminosity Lp =16 
mJy kpc^ labelled by survey number as given in Table |I|. (Bottom:) Maximum survey distance, 
Dmaxi averaged over direction for each of the eight surveys and for a luminosity of 16 mJy kpc“^. 
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The ratio of V to Vmax from Eq. may then be written as 

E 

^max 

Applicaton of Eq. ^ involves subtleties that depend on whether the flux density reported for 
a given object is influenced by the random process associated with DISS. DISS generally increases 
the volume in which a pulsar can be detected (cf. Appendix]^). Saturated DISS that is not 
quenched by time-bandwidth averaging modulates the flux density by a random variable drawn 
from an exponential probability density function (pdf). Though the modulation is less than unity 
more often than not, the net effect is to increase the volume by a factor r(5/2) ~ 1.33. 

Figure ^ shows E/Enax for field MSPs plotted against \z\ = Zlsin|6|. Horizontal error bars 
reflect uncertainties in the measured flux density and distance errors, the latter also determining 
the vertical error bars on \z\. In most cases, we have used the 400 MHz flux density reported 
by Taylor, Manchester &: Lyne (1993), which is usually an average of many observations and is 
influenced minimally by DISS. Flux calibrations are typically only about 20% accurate. Weak 
pulsars can appear brighter than average due to DISS at the time of discovery, however, so 
the discovery flux > Smin while S < Smin- For these cases (J0034-0534 and J0711-6830), 
E/Enax ^ 1- The intrinsically brightest MSP, B1937+21, is detectable to only ~ 8 kpc despite its 
large luminosity Lp>3100 mJy kpc^ because, at its low galactic latitude, dispersion and scattering 
effects grow rapidly with distance. Consequently, we argue that the claimed upper distance limit, 
Du 15.7 kpc, is a factor of two too large. The moderate-latitude pulsar J1643-1224 is attributed 
only a lower bound on its distance by the TC model, D > 4.8 kpc, because its DM cannot be 
accounted for by the model. We suspect that this pulsar’s DM is enhanced by unmodeled ionized 
gas along the line of sight and that the distance is most likely less than 4.8 kpc. In the absence of 
further data, however, we use the distance lower bound as is. 

Apart from the pulsar with a questionable distance estimate (J1643-1224), Figure ^ shows 
that MSPs are to be found at only low values of \z\, suggesting, therefore, that the scale height for 
MSPs is ~ 0.5-1 kpc. To properly estimate the scale height requires careful accounting of selection 
effects in MSP surveys, as we do in §Q. However, Arecibo surveys at high latitudes search to 
several kpc for typical luminosities. The absence of high \z\ pulsars is therefore especially striking. 
The Arecibo MSPs also tend to have small values of VjVmax, as would be expected for surveys 
that search well beyond the scale height of the population. 




S 


3/2 


(5) 


4. LIKELIHOOD ANALYSIS 

4.1. Observables, Assumptions and Statistical Method 

A survey for MSPs typically searches many beam areas for each MSP discovery. The spatial 
distribution of MSPs determines this yield, along with the survey sensitivity as a function of 
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Fig. 4.— A plot of \z\ vs. VfVmax for field millisecond pulsars. Filled circles denote pulsars 
discovered at Arecibo. Filled squares indicate MSPs found at Parkes and Jodrell Bank. The filled 
triangle denotes the lower bound on z and VfVmax for J1643-1224, whose distance estimate is only 
a lower bound, as discussed in the text. 
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period, the period distribution, and the luminosity function. Here we derive the likelihood 
function for a survey, taking these factors into account. We take as observables the directions 
of all beam areas searched, the survey sensitivities in these directions, and the parameters that 
describe individual pulsars, including direction, period, flux density, and dispersion measure. We 
also use distance estimates based primarily on the electron density model of TC. Such distances 
are imprecise and our method takes into account the large uncertainties in distance that translate 
into large uncertainties in implied luminosity. 

Consider a telescope beam with solid angle fi/,. The mean number of MSPs expected in the 
beam per unit period, luminosity and distance is 

= ni.0""p(0.(6) 

where the number density of MSPs, np{D, i,b), is an arbitrary function of position. We have 
assumed that the joint probability distribution of period P, luminosity Lp and position is 
factorable. The physical assumption is that the distribution in space is independent of the 
distribution of intrinsic pulsar properties {P and Lp) and, furthermore, that P and Lp are 
uncorrelated. We write the period pdf, fp{P), and luminosity pdf, /ppi^p), each with unit 
normalization. With the small number of MSPs currently available there is scant evidence that 
the factorization is or is not appropriate; however, in the future, our method can be applied easily 
to more complicated joint distributions if warranted. In particular, we defer to another paper 
exploration of the joint statistics of Lp and P. Below, we specialize to disk and disk + diffuse 
models and we take into account the variation of the telescope gain across the its beam. 

To calculate the mean number of MSPs expected per beam, we integrate Eq. to obtain 

/P C TTLCbX 

dPfp{P) J dLphpiLp) dDD^np{DJ,b). (7) 

The volume searched per beam, averaged over P and Lp, is 

SVs = ^l dPfpiP) I dLpfL^{Lp)Dl,,. (8) 

Survey sensitivities are implicit in Dmax, as discussed in §|^ For detections, we take into account 
the constraints that exist, due to post-discovery observations, on period, flux density, and distance: 
P ± AP/2, S ± AS/2, and D G [Dl,Du]. Integrating over the subvolume bounded by these 
constraints, the mean number of MSPs is 


{^p)d 


[ dPfp{P) 

Jp±AP/2 


Qb f dPfp{P) 

Jp±AP/2 


(•min[D[/,Drj 


iDr 


dDD%{D,£,b) I dS'fppiS'D^'' 

Js±AS /2 


( 9 ) 


rS^Dl 

IsiDl 


dLph^iLp) 


rram[Du,Dmaa:,iLp/Siy^'^] 

lma.x[DL,{Lp/S2yt2] 


dDD‘^np{D,i, b). 
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where Si ^2 = <5 =F AS'/2. For each known pulsar, D^ax is the maximum distance for the survey 
that could have detected each pulsar. Several MSPs in our analysis were first discovered in other 
surveys but were subsequently detected (rediscovered) in one of the eight surveys and we analyze 
them accordingly. For a given beam, the Poisson probabilities for detecting zero or one MSP are 

Pi = (iVp)e-W. (10) 


We construct the survey likelihood function as the product of nondetection (ND) and detection 
(D) factors; 

C- = Cnd^d, ( 11 ) 

where, for N^, total beams searched, Mp MSPs found, and assuming (Np) <C 1, 


Cnd = exp [-X] (IVp). 1 

(12) 

\ i=i / 


Mp 


C-'D = {^p) Dk 

k=l 

(13) 


The log likelihood is 

Ni Mp 

A^inC = -Y, {Np)^ + X: (14) 

j=l k=l 


The likelihood function may be simplified if we factor the pulsar number density into a 
constant no times a shape factor: 

np{D,i,b) = noh{D,i,b), (15) 

where h{D,i, b) is dimensionless and has a maximum of unity. Substituting, the likelihood function 
becomes 


Mr, 


A{0, no) = Mpinno — noVd + X] 


(16) 


k=l 


where the vector 0 denotes the set of parameters other than no- We define the survey detection 
volume as the sum over beams, 

Nb 

Vd = J25Vdj, (17) 

i=i 

where (dropping beam labels) 

d{Np) 


SVd = 


dno 


and the constrained subvolume per discovered MSP is 

sv, = 

dno 


(18) 


(19) 
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The survey detection volume Vd is the volume searched weighted by the dimensionless MSP 
space density, h. The expected number of MSPs in a survey is simply noV^. Eq. E applies to a 
single-component density model, such as the disk distribution we consider in the next two sections. 
Multiple components require the alternative treatment of §^. 


Maximizing C with respect to no, we obtain the best fit number density (for a specific set of 
parameters, 6) 

ho = (20) 


Vd 


Substituting, the log likelihood becomes 


Mp 

A(0, ho) = Mpin ho — Mp + ^ in 5Vp. 

k=l 


( 21 ) 


For no 7 ^ ho, the variation in the log likelihood is 


A(0, no) - A(0, ho) = Mp 


ini"^]- 
no 


no - no 
no 


Mp / no - ho 

2 V ho 


( 22 ) 


where the approximate, quadratic form holds for |no — ho|/ho 1 . 


We want to know the marginal distribution of each parameter. For a given parameter 6j G 6, 
the marginal pdf is the normalized integral over all other parameters 


feM) = 


lexc-e, I dnoC{e, no) deC{e, ho)ho 


f dO JdnoC{9,no) 


f d6C{6, ho)ho 


(23) 


where the integral subscript ‘exc. Oj’’ means that all parameters except the one are integrated 
over. The approximate form in Eq. 23 assumes a sharp peak about ho and becomes an increasingly 
good approximation as Mp grows. The marginal pdf for no is 


fnoino) — 


JdeC(e,no) 

J dO JdnoC{6,no) 


Jde£{9,ho)e 


Mp f no - no Y 


no 


f d9C{9, ho)ho 


(24) 


For disk models the areal, or column, density of MSPs is less model dependent than the 
number density and scale height separately. The column density is A^o = dno^^z, where 7 / is a 
dimensionless factor of order unity and is a scale-height parameter. The pdf of No is calculated 
by marginalizing £ over all parameters except no and az- The resultant joint pdf /no.o-^ is then 
integrated according to: 

fNo{No)=V~^J daz<T~^fno,aA^o/Wz,(^z)- (25) 

For disk models considered below, rj = 2 for an exponential in z and rj = for a Gaussian in z. 
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4.2. Telescope Gain 


MSP surveys usually involve drift scans or sustained pointings toward specific sky positions. 
The telescope’s gain toward a given source varies over the analyzed portion of the drift scan and 


is a function of the source’s position relative to the beam center (see, e.g., Camilo, Nice &: Taylor 
1996 ). We account for gain variations by replacing the beam solid angle in Eq. with a sum 
over equal solid-angle terms 


Qh 




(26) 


m=l 


where the telescope gain varies with m, Gm- The minimum detectable flux density Smin is 
therefore a function of m. For some drift-scan surveys, we take into account that the data are 
analyzed in data blocks that overlap by some fraction (usually 50%). 


For drift scans, G varies with time over the data set and the offset from the beam center in 
declination is also taken into account. The sum in Eq. |2^ becomes a sum over discrete steps in 
declination. For pointed (tracking) observations, we use actual pointing directions and break the 
beam into equal-solid angle annuli about the beam center, which we sum over as in Eq. |^. We 
find that only a small number of subbeam elements is needed to account for the shape of the 
beam, e.g. ~ 2 or 3. 


4.3. Interstellar Scintillations 


In Appendix ^ we derive the effects of diffractive interstellar scintillations on flux densities 
and on the (pseudo) luminosity function. To use these results, we replace in Eq. ^ with 
the corresponding ‘scintillated’ luminosity function, as defined in the Appendix. We do so 
for surveys assuming that specific sky positions are observed only once. However, we use the 
unscintillated luminosity function in Eq. 10 because flux densities reported for the known pulsars 
are generally long-term averages of many independent measurements. 


4.4. Comparison with Other Statistical Methods 

Our statistical method differs substantially from other studies of the MSP population. A 
common approach to population studies, including pulsars and gamma-ray bursts, makes use of 
nonparametric estimators. The rationale is to try to draw inferences about certain properties of 
the population without assuming a specific class of models. In contrast, our likelihood analysis 
makes very specific assumptions about the class of models to be examined; for example, we 
have assumed a priori that all probability distributions are continuous. The differences between 
parametric and nonparametric treatments highlights some of the strengths and weaknesses of our 
approach. 
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It has been shown (Loredo and Wasserman 1995) that nonparametric estimators may be 
derived from a special maximum likelihood model solution. Since our parametric treatment is also 
based on a maximum likelihood analysis, it is straightforward to study the relationship between 
alternative methodologies by making a comparison of the assumptions made in the two searches. 
The special solution leading to the nonparameteric estimators of interest comes from a search 
for a maximum likelihood solution amongst all functions and generalized distributions (i.e. delta 
functions) with equal a priori weight. This class of functions is so large that the most likely 
model is always one which exactly and precisely describes the observed data; thus, nonparametric 
estimators satisfy the rationale for which they are introduced. This contrasts with the parametric 
treatment adopted here for which the class of functions is (by comparison) extremely small. We 
liken nonparameteric estimators to models with large numbers of free parameters. 

In deciding what treatment to adopt, it is helpful to appeal to the Bayesian odds ratio to 
decide whether adding a new parameter to a model is justified by the better description of the 
data it may entail. Roughly speaking each newly added parameter will improve the quality of the 
model’s description of the data. The odds ratio allows a quantitative decision to be made whether 
to adopt the more complex model by weighing the improvement in the description against the 
additional freedom to fit arbitrary data sets. The situation for the MSP’s is that the population 
is rather small and we have anticipated (without any detailed investigation) that the odds ratio 
will favor models with relatively small numbers of parameters. We have therefore focused in this 
paper on parameteric methods with small numbers of parameters. 

An additional factor in our choice of parametric methods is that it is straightforward to 
include ancillary information about the population (e.g. continuity of the model), whereas in 
nonparametric approaches such constraints are difficult to incorporate. Moreover, we find the 
parametric approach naturally allows the inference of population parameters of significant interest 
(e.g. cutoffs in the period distribution). 

The main drawback of the parametric approach is that the results apply only to the particular 
set of models that the parameters can describe. If the real data were much better described 
by some completely different unstudied model, one would have no indication of that fact. In 
this paper we have considered several plausible models but these cannot begin to describe all 
possibilities. 


A number of pulsar population studies are based in whole or in part on such estimators 
(Vivekenand Sz Narayan 1981, hereafter VN; Phinney Sz Blandford 1981 and Narayan 1987| ). To 
be a bit more descriptive, in the VN method a scale factor is calculated for each object detected 
in a survey. The factor represents how many pulsars with the same period P and luminosity Lp 
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exist in the Galaxy given the fraction of the Galaxy searched. In our notation, the scale factor is 


S{P,Lp) 


roo 

^ Ob, / dDD^np{DJ,b) 

full sky 


Ni 


i=i 



dDD^np{D,i,b) 


(27) 


where D^ax, as before, depends on many survey and pulsar parameters, including P and Lp. The 
number of pulsars in the Galaxy is then calculated through a sum over detected pulsars as 


Ngal- E 
i=l 


(28) 


The resultant total number of pulsars is a mean value similar in nature to the mean value of the 
number density, ho, that we have calculated. One drawback is that the VN method estimates the 
number of pulsars in the Galaxy exactly like those actually detected. In other words, it explicitly 
includes contributions to the mean only at the periods and luminosities of the known pulsars. It 
is inherently discrete as compared to our likelihood method based on continuous distributions. 
Another drawback is that the method does not directly allow computation of confidence intervals. 
Finally, since the scale factors are calculated only for the detected pulsars, there is no means for 
estimating the cutoffs of the distributions of P and Lp. Below we compare our results on MSPs 
to those of Lorimer et al. (1995) and Bailes & Lorimer (1995; hereafter BL) with these issues in 
mind. 


5. DISK MODEL FOR MILLISECOND PULSARS 

5.1. Method 

The simplest spatial model is a disk with constant scale height so that the density is a 
function of z only. We let np{z) = nc[hp{z) with hp{0) = 1, where Ud is the midplane density that 
corresponds to no in §^. 

The parameters to be solved for describe the period, luminosity and z distributions, fp{P), 
fipiLp) and np{z). We have considered three models for hp{z): (1) a Gaussian function in z 
with an rms value of z given by (2) an exponential model with 1/e scale height and (3) a 
numerically derived distribution of NS orbits, neither Gaussian nor exponential in form, discussed 
in §^. For the luminosity and period pdfs, we adopt power-law functions, i.e. fpp oc Lp and 
fp oc with respective lower and upper cutoffs, Lp,,Lp 2 and Pi, P 2 . 

The greatest computational effort goes into calculation oi Cpip (Eq. |T^). We computed it 
efficiently by summing the D integral in Eq. over the survey beam areas for a grid of cr^, P, and 
Lp] as stated before, we use the scintillation-modified luminosity function in this computation. 
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Next we form the likelihood for different model parameters aLp,ap, Lp^, Lp^, Pi, P 2 by calculating 
integrals over P and Lp with weights fp and fpp (cf. Eq. |^. 

We maximized A by varying the parameters (or subsets of the parameters) over a grid. 

We kept the upper cutoffs on the period and luminosity distributions fixed at P 2 = 20 ms and 
Lp 2 = 16,000 mJy kpc^. The period cutoff corresponds to the selection used to define the sample. 
Since the number of objects decreases rapidly as P increases, the upper cutoff plays little role in 
any of the results below. The luminosity cutoff corresponds to the maximum possible luminosity 
in the observed sample. We also tested the effects of varying Lp^ and found that results are not 
sensitive to this parameter. Exclusion of B1937+21, the most luminous MSP, allows a much 
smaller value for Lp^ to describe the remaining 21 pulsars in the sample; but none of the other 
results below are substantially altered. 


5.2. Results 

The five parameters (Pi, Lp^, aLp,ap, and a^) were varied over a grid to find the maximum 
A. We formed marginal pdfs according to Eqs. and 
Eigure shows the marginal pdfs for each of the six parameters (the above-mentioned five and the 
number density, n^). Using these pdfs, we calculated the confidence intervals on the parameters 
that are given in Table The maxima are well-defined and easily located. 


24. Results are summarized in Table 


5.2.1. Minimum Period Pi and Period Distribution Slope ap 

Naturally Pi must be less than or equal to the period of the shortest-period MSP in our 
sample. When other parameters are held fixed, it is straight forward to show that A must 
decrease as Pi is made smaller. The best-fit, minimum period lies only slightly below that of the 
most-rapidly-spinning, known pulsar, B1937-I-21 (1.56 ms). However, the data allow Pi < 1.56 
ms at a reduced level of confidence. The results are given in Table ^ for both the Gaussian and 
exponential models. The cutoff is > 1 ms at 95% confidence and > 0.65 ms at 99% confidence. 

The period distribution falls off steeply with period, implying the existence of many objects 
at small P {dN/dP oc /p oc p-2.o±o.33^^ jg known that physical instabilities will act on 
neutron stars with very short rotation periods. Ignoring the magnetic field and assuming accretion 
from an inner edge of a Keplerian disk. Cook et al. (1994a,b) have shown that 1.4M0 neutron 
stars can be spun up to critical rotation periods (well under 1 ms) for a variety of equations of 
state without triggering radial instability, e.g. exceeding the maximum neutron star mass. (The 
results do not assure stability against non-radial modes and the associated gravitational wave 
emission.) Our overall fit for the period distribution suggests the existence of MSPs faster than 
those that are currently known (1.56 ms) in view of the fact that the theoretical stability analyses 
do not rule out such objects. Of course, there may be evolutionary reasons that such objects do 



logic fp, logjo f logjo f 


- 19 - 




Fig. 5.— Marginal pdfs for the 6 parameters of the exponential disk model for the millisecond 
population. is the column density (kpc“^) of MSPs, which we show instead of the number 
density, n^. 
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not occur and we discuss the significance of the cutoff Pi next. 

The specihc value of Pi depends, of course, on our assumption of a power-law distribution for 
P. We have not explored other mathematical forms, but reasonable alternatives include a power 
law that flattens for periods less than some critical period and cuts off at Pi or a distribution that 
rises slowly from zero at Pi and peaks at or near 1.56 ms and then follows a power-law form like 
that we have fitted. It is easy to see that such alternative period distributions will lead to smaller 
Pi than we have derived. The reason is that they imply that smaller volume has been surveyed 
for P < 1.56 ms, so the allowed Pi can be smaller. Therefore, our derived Pi using the power-law 
distribution is a maximally allowed value and suggests, conservatively, that the period range for 
MSPs may extend to as small a value as 1 ms (95% confidence) or 0.65 ms (99% confidence). 

Harding (1984) analyzed the slope of fp, assuming a steady-state flux with births balanced by 
pulsars crossing the Hubble line. She showed that if pulsars are born with a powerlaw distribution 
of B (oc B^) and with initial period P approximately oc B (accretion spin up models imply P®/”^) 
then the resultant fp oc P^. Today it is known that the spin-down times for the observed MSPs 
are too long for a steady-state to be attained. However, with similar assumptions we find the same 
slope in the period range [Pmin{Th), Pmax], where Pmin(Th) is the period reached after a Hubble 
time (Th) by the minimum initial period object (minimum magnetic field) and Pmax is the longest 
period at birth. (Different slopes are found in other period subintervals. Additional discussion will 
be found in Chernoff & Cordes 1996a). Thus, one possible interpretation of the steep period slope 
is that the field distribution oc iJ-2.o±o.33^ However, it is difficult to derive robust constraints on 
the field distribution without knowing both the time dependence of magnetic fields during the 
spinup process and the spindown law for MSPs subsequent to the spinup phase. 


5. 2.2. Scale Height 

The inferred Gaussian scale (0.65 kpc) and exponential scale (0.50 kpc) are in rough 
agreement. The values indicate that the MSPs have a relatively small scale height, comparable to 
the oldest disk stars. Though the confidence intervals overlap, the actual shape of the distribution 
plays some role in the value of the scale height parameter and motivates, in part, a more physical 
analysis based on motion of objects in the Galactic potential (§^). 


5.2.3. Minimum Luminosity Lpi and Slope app 

The luminosity pdf of our best fit, dN/dLp oc /ppiLp) oc is similar to that of 

long-period pulsars (e.g. Lyne, Manchester & Taylor 1985). Total numbers are dominated by 
weak sources. The lower cutoff is Lpi = I.I^q^ mJy kpc“^ and is largely determined by the 
absence of nearby sources. We have shown for long-period pulsars that fpp is strongly influenced 
by geometrical beaming effects, the distribution of true luminosities, the spin down law and a 
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death line (Chernoff & Cordes 1996a). Because all four of these elements may differ between 
high-held pulsars and MSPs, we currently regard the similarity between the long-period and MSP 
luminosity distributions as fortuitous. 

In the past, the disk-determined fipiLp) (slope and cutoff) has also been used to make 
inferences about the number of MSPs in globular clusters. On evolutionary grounds, many 
properties of disk and globular cluster MSPs might be expected to differ (e.g. distributions of 
luminosity, spin period, orbital period and velocity). Since the nearest cluster is too distant to 
allow direct measurement of the luminosity function near Lpi, usage of the disk-determined form 
is necessary for many purposes. Fruchter &: Goss (1990) measured the radio ffux from nearby 
globular clusters and estimated ~ 10^ MSPs in the Galaxy’s globular cluster system. Our best 
ht luminosity distribution, with cutoffs, is consistent with the one they assumed and does not 
alter the size of this estimate. Likewise, estimates by Foster Sz Tavani (1992) and Johnston, 
Kulkarni &: Phinney (1992) of the shape of the luminosity pdf for MSPs in globular clusters are 
also consistent with our best-fit for disk MSPs, though both groups were unable to determine 
the lower luminosity cutoff and, hence, the absolute normalization. Wijers & van Paradijs (1991) 
find far fewer globular cluster MSPs than do Fruchter &: Goss or Johnston, Kulkarni & Phinney 
even though they adopted a lower luminosity cutoff three times smaller than that of Fruchter 
&: Goss; the difference is probably related to their assumed dependence of luminosity on spin 
period and spin period derivative that was based on young, high-field pulsars. Analysis of globular 
cluster MSP populations should probably use a treatment similar to this paper’s but applied to 
cluster-only data. 


5.2.4- Correlations 

Most of the derived parameters are uncorrelated. However, ap and Pi are positively 
correlated as are and Lp.^, while is negatively correlated with the lower cutoffs in period 
(Pi) and luminosity (Lpj^). Figure^ shows contours of constant likelihood plotted against pairs of 
parameters while holding all other parameters fixed at values that yield the maximum likelihood. 


6. DYNAMICAL MODELS 
6.1. Birth Kick Determination 

In §11 we assumed functional forms for the 2 : distribution of MSPs and fit for the associated 
scale height parameters. These parameters describe the present-day MSP distribution without 
regard to the orbit about the Galaxy. We have constructed a dynamical model that connects 
“birth parameters” to today’s spatial distribution as follows. We model the birthrate density of 


Lp^ (mJy kpc o(p 



0.5 1 1.5 0.5 1 1.5 2 

P, (ms) (kpc) 


Fig. 6.— Selected contour plots of the log likelihood for the exponential model plotted against 
pairs of parameters while holding the other four parameters fixed at their values that yield the 
maximum likelihood. Contour spacings are unity in natural log units and the first contour is a 
factor 1/e from the peak. 
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MSPs 


h{R, z) = g{R) exp 



(29) 


where R is (cylindrical) Galactocentric radius, z is height about the plane, az^b is a scale height 
parameter and g{R) is a surface density function, taken to be either constant (“uniform model”) 
or exponential with scale length 3.5 kpc (“exponential model”). 


The birth velocity is the circular rotation velocity plus a kick component. Note our use of 
“kick” includes any momentum impulse imparted to the pulsar’s progenitor or companion, if 
in a binary. The angular distribution of the kick is isotropic and the velocity magnitude has a 
distribution oc . After birth, the MSP trajectory is determined by integration of the 

orbit about the Galaxy in a simplihed model of the gravitational potential (Pacyzynski 1990). We 
ignore the role of scattering from irregularities (e.g. GMG’s, spiral density waves, massive black 
holes) in the calculated motion. We first discuss the kinematic properties of the MSP population 
inferred from the smooth model and next assess the degree to which our conclusions may be 
modified by the diffusion of stellar orbits. 

About 4 million orbits were integrated over time spans of 10® years, sufficiently long that the 
derived vertical distribution was stationary and well-mixed. For specific birth parameters, the 
vertical distribution of MSPs in the vicinity of the Sun (e.g. in an annulus of Galactocentric radii 
from 7.5 — 9.5 kpc) was calculated by appropriately weighting and combining the results for the 
individual orbits. 


The statistical analysis described in previous sections was carried out to determine the birth 
parameters {ay and the intrinsic pulsar population parameters) of the uniform model. The results 
(Table ^ give a peak value ay = 52^(( km s“^. The initial scale height is not well-determined 
by the data and was held fixed, az,b = 0.1 kpc. The column density of MSPs was calculated 
from Eq. |^. The result is consistent with values obtained using the Gaussian and exponential 
spatial models. Figure ^ illustrates the density distribution vs. z, comparing the range of allowed 
exponential fits (§|^ to the most likely dynamical model. The differences are subtle and suggest 
that the assumed exponential form should be an adequate local description for many purposes. 


6.2. Kinematics of Today’s Population 

Kinematic properties of the MSP population may be inferred from the dynamical model. 
For example, the distribution of parallel and perpendicular velocities relative to the LSR are 
easily derived from the orbital calculations. Figure ^ shows the distributions for all simulated 
objects within 1 kpc of the Sun for the most likely dynamical model. The expected transverse 
motions are small; approximately 99% of the MSPs have AV± < 150 km s“^. As MSP samples 
increase in number, detailed distributions like these will provide important additional constraints 
on modeling. 
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Fig. 7.— Comparison of the z distributions for the best fit velocity model (noisy solid line) and 
three exponential models with scale heights of 0.41, 0.55 and 0.78 kpc (dashed lines). 
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Fig. 8.— Histograms of observed transverse speeds (solid) and line-of-sight (dashed) velocity for 
the best fit velocity model of §^. 
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Next, we consider the velocity ellipsoid of MSPs. Let vr, vt and Vz be the components of 
velocity in the cylindrical radial direction, in the tangential direction (parallel to the local circular 
velocity, e.g. I = 270° at the solar position) and out of the plane, respectively. In the well-mixed 
state, the most interesting non-zero moments are < v\ >, < >, < vl > and < vt >■ Table 

^ lists the first and second moments for a range of MSP birth models with ay = 20, 40, 60, 

80 and 100 km s“^, for two scale heights az^h = 0.05 and 0.15 kpc, and for the uniform and 
exponential surface density distributions. (All velocity moments are given in units of ctv-) When 
the kick velocity is small compared to the rotation velocity and when disk properties do not vary 
significantly over the range of radii sampled, the results of epicyclic theory are directly applicable. 
For local objects < > / < v\ >= \B/{A — B)\ has the observed value 0.45 ± 0.09 (for Oort 

constants A = 14.5 ± 1.5 km s“^ kpc“^ and B = —12 ± 3 km s“^ kpc“^ 

1987[ ). The model-calculated value of 0.45 at ay = 20 km s“^ is in good agreement with the value 
inferred from the observed Oort constants. Here, we will concentrate on the changes that occur as 
ay increases and that are indicative of some of the differences between the velocity distributions 
of MSPs and of disk stars. A global model is necessary since a local epicyclic treatment for the 
MSPs is not well-founded. For example, with a kick of 60 km s“^ particles observed at the local 
position could come from initial radii in the approximate range (0.5 — 2.6) x Rq, where Rq is the 
Sun’s distance from the Galactic Center. Also, kicks of this size create vertical excursions of > 0.3 
kpc, spoiling a harmonic approximation to the potential. 

Table shows how the basic moments change as ay increases. We briefly note the most 
important conclusions: (1) A clearly noticeable effect is the occurrence of a non-zero tangential 
motion measured with respect to the local circular velocity (“asymmetric drift”). For ay = 60 km 
s“^ the magnitude is ~ 13 km s“^ in the uniform model (~ 25 km s“^ in the exponential model), 
an effect that is potentially detectable in a relatively small sample of objects with well-determined 
velocities. (2) The velocity ellipsoid (with axial ratios \J< >:\J< uf >.\/< >) becomes 

rounder as the magnitude of the kick grows. (3) The birth distribution in Galactocentric radius 
affects the value of all the non-zero velocity moments including the shape of the velocity ellipsoid 
and the magnitude of the asymmetric drift. (4) The imprint of the birth scale height is essentially 
absent for objects with ay > 60 km s“^. 

Determinations of asymmetric drift would provide valuable information on the birth locations 
of MSPs. Proper estimation of the effect will require more field MSPs than are currently known 
and careful treatment of distance errors. We defer a detailed discussion to another paper. 


Binney &: Tremaine 


6.3. Orbital Diffusion 

The model calculations presented above assume a regular background potential. Older stars 
are well known to have larger velocity dispersions, presumably from interaction with small-scale 
fluctuations in the gravitational field, but the actual physical source of the irregular field is not 
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well understood ([Wielen 1977]) . The oldest stars, K and M giants of age 9 x 10® yrs, reach total 
dispersions of 77 km s“^; for comparison, using interpolated values for the uniform model in Table 
^ we estimate that the best fit model for the MSPs (cry = 53 km s“^) implies a total dispersion of 
84 km s“^. In fact, the MSPs suffer comparable energy input from kicks and from diffusion. The 
key assumption is that the MSP population includes members with ages ranging uniformly up to 
the age of the Galaxy, so that the average effect of diffusion will be less than it is for the oldest 
stars. Using the velocity dispersion data of K and M giants with ages (0.3 — 9) x 10® ( Wielen 
1977 ), averaging uniformly in time, we infer that the root mean square dispersion is ~ 50 km s“^. 
We suggest that the residual dispersion of 67.5 km s“^ (i.e. VSd® — 50^) is due to kick(s) unique 
to MSP evolution. This 3D dispersion would then correspond to a ID kick of ~ 39 km s“^. 


6.4. Conclusions 

The best fit uniform model implies ay = 53 km s“^; this is an upper limit because 
gravitational scattering processes are ignored in its estimate; the scale of the kick is ~ 40 km s“^ 
assuming MSPs are long-lived and born at a uniform rate. If MSPs are visible for less than a 
Hubble time, the kick size will increase; if most of today’s MSPs were formed early in the Galaxy’s 
life, the kick size will decrease. 


7. DISK -h DIFFUSE MODEL 
7.1. Method 

The MSP distribution may be more complex than a single disk component with small scale 
height. For example, there may exist a population of MSPs that fill a halo-like region around the 
disk. 


If MSPs are distributed in two components, the log likelihood becomes 


A{e,nd,nh) 


Mp 

-[ndVd + nhVh] + ^in 

k=l 




(30) 


Here, we label disk quantities with ‘d’ while ‘h’ denotes diffuse (halo-like) contributions; we 
suppress the dependences of the volumes on other parameters. Maximizing A with respect to 
and Uh, we hnd that the best-fit number densities hd and h/j satisfy 


ndVd + fihVh = Mp. 


(31) 


Also, if we take the = 0 case as a hducial solution, which is our result in for the disk-only 









- 28 - 


model, the log likelihood for Uh ^ 0 may be expanded as 


Mp 

A{e, Ud, Uh) = A(0, Ud, 0) - rihVh + 

k=l 


V udSvi^^) ■ 


(32) 


7.2. Diffuse Models 


One might expect a diffuse distribution of MSPs for any of several reasons. (1) The probability 
distribution of birth velocities may extend to values much larger than typically allowed by the 
assumed Gaussian or exponential forms. The high velocity MSPs would oscillate to higher z 
distances or escape the Galaxy all together. (2) MSPs born in globular clusters may be ejected by 
dynamical interactions or when a cluster is tidally dissolved. Such objects would have a spatial 
distribution like the parent systems assuming the ejection velocities were small compared to the 
rotation velocity. (3) Spheroid stars may evolve and produce long-lived MSPs just like disk stars 
(e.g. by accretion-induced spinup). Such objects would have a spatial distribution like the Pop 
II spheroid. (4) If the formation of the Galaxy involved hierarchical merging of smaller objects 
containing disk-like structures, their MSPs will be cannibalized. Such objects might follow the 
dark matter halo distribution. 


Without further considering the merits of these basic scenarios, we will adopt several 
geometrical distributions for the the putative diffuse population and place upper limits on the 
number densities. Consider a density model for MSPs of the form 


nh{r) 


rih 


1 + (r/rh)'^ 


-Sh/2 


(33) 


where r is the radius from the center, rh is the characteristic radius, and Sh is the power-law index. 
Taking s/j = 0 gives a uniform density halo, our reference model (in practice all distributions are 
truncated at 50 kpc). Taking s/j = 2 and rh = 5 kpc gives an isothermal distribution with large 
core. Taking Sh = 3.5 and r/j = 1 kpc gives the observed globular cluster distribution (Thoma^ 


19891) . 


Using these models, we calculate the diffuse pulsar density by integrating Eq. and we 
evaluate the halo volume factors Vh and 6Vp^ (Eqs. ^ and [l^ ). We combine these with the 
analogous disk quantities and examine a grid in Ud and Uh to find the distribution of likelihood 
values. 


7.3. Results 

Eor our reference model, we find that the pure disk model is favored by a huge factor implying 
an upper bound on the diffuse density from the fitting is Uh ^ 0.4 kpc“^ (90% confidence, cf. 
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Table D- Figure ^ shows likelihood contours for the disk and diffuse densities, with a maximum 
at rifi = 0. The marginalized densities are shown in Figure 0 for the uniform density model. 

For the other two models, we have expressed the results in terms of limits on the density 
parameter re/j (the value at the center of the Galaxy) and, equivalently, on nh{Ro) where Rq is 
the Sun’s galactocentric radius. Though we have made calculations explicitly for the nonuniform 
density models, the local values are close to the reference model values. This follows because on 
the Galactic scale most surveys probe regions near the solar system. 


7.4. Disk, Spheroid, Halo and Globular Cluster Contributions 


The local column density of MSPs, ~ kpc may be combined with the disk 

surface mass density (~ 66 ± 8Mq pc“^ for Oort K giants, Bahcall 1984|) to infer that the number 


of MSPs per unit disk mass {dN/dM)disk ~ 7 . 6 ^ 3 ;! x (the range reflects only the 

uncertainty in N^). The total number of MSPs in the Galactic disk scaled to the disk mass M^isk 
is S.Ol];;® X l{)‘^{Mdisk/^ x lO^^M©) (for example, M^isk = 3.7 x W^^Mq |Bahcall &: Soneira 1982 ] 
by one estimate; (3.5 — 4.6) x [paldwell fc Qstriker 198l| by another). The total does not 

include a correction for beaming. 

Estimates of the local spheroid mass density are uncertain, e.g. psph = 1-88 x lO“^M 0 pc“^ 
( Bahcall, Schmidt &: Soneira 1982| ) or psph = (1.11 — 1.25) x lO“^M 0 pc“^ ( Paldwell &: Qstriker 
1981). If {dN/dM)disk = {dN/dM)sph, then the spheroid makes a contribution to the MSP 
number density Usph = 0.14 kpc“^ or (0.84 — 0.95) kpc“^, respectively. The upper limit we have 
derived for a uniform density model, <0.42 kpc“^, is marginally consistent. Future observations 
should be able to constrain contributions to the MSP population from Population II progenitors 
more strongly. 


Estimates of the dark matter halo density are phaio = 9 x 10 ‘^Mq pc ^ ( Bahcall, Schmidt &: 
Soneira 1982 ) or phaio = (5.9 —10.2) x 1O“^M0 pc“^ ([Caldwell h Qstriker 198l]) . If the dark matter 


halo component satisfied {dN/dM)disk = {dN/dM)haio then its contribution is Uhaio = 6.8 kpc“^ 
or (4.5 — 7.7) kpc“^, respectively. Qur limit on implies {dN/dM)halo/{dN/dM)disk < 0.06 or 
(0.05 — 0.09), respectively. 


Today’s globular clusters are known to have a significant enhancement of MSPs 
relative to the disk. With considerable uncertainty, Phinney and Kulkarni (1994) estimate 
{dN/dM)gc ~ 50{dN/dM)disk- If half of the original globular cluster system has been destroyed 
(e.g. a total mass M ^ hx Mq), if the MSP content was similarly enhanced and if these MSPs 
orbit like the observed clusters, then the contribution to the local mass density is 1.1 x 1 O“^M 0 
kpc“^ and the MSP number density is 4.2 x 10“^ kpc“^. Qur limit on does not provide a 
strong constraint. 
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7.5. Conclusions 

The observations place an upper limit on a diffuse halo-like contribution to the MSP density 
that is roughly 1% of the MSP disk density at midplane. 


8. SPACE VELOCITIES OF MSPS 

Our results indicate that millisecond pulsars are a low-velocity population, at least when 
compared with young, high-held pulsars. We have found that the 3D rms velocity of MSPs in the 
galactic disk is ~ 84 km s“^ , about a factor of 5-7 lower than that of young, strong held pulsars 
(Lyne & Lorimer 1994; Cordes &: Chernoff 1996). We have reached this conclusion by determining 
the spatial distribution of MSPs, by excluding the existence of a signihcant non-disk population 
and by modeling the motion of objects in the gravitational potential of the Galaxy. 


8.1. Comparison with Proper Motion Data 

We may compare our results with direct measurements of proper motion using interferometric 
and pulse-timing methods; the indirect method of interstellar scintillation has also yielded 
determinations of MSP transverse speeds. To date, there are timing proper motions on eight 
MSPs: J0437-4715 (Bell et al. 1995), B1257-I-12 (Wolszczan 1994), J1713-I-0747 (Camilo, Foster &: 
Wolszczan 1994), B1855-I-09 &: B1937-I-21 (Kaspi et al. 1994), B1957-I-20 (Arzoumanian, Fruchter 
&: Taylor 1994), J2019-I-2425 and J2322-I-2057 (Nice &: Taylor 1995). There are also scintillation 
speeds on some of these and other pulsars, B1855-I-09 (Dewey et al. 1988), B1937-I-21 (Cordes 
et al. 1990), and J0437-4715, J1455-3330, J1730-2304 & 2145-0750 (Nicastro & Johnston 1995). 
These MSPs have transverse speeds that are less than 100 km s“^ , except for B1257-I-12, which 
has a speed of 285 km s“^ at its nominal distance of 0.62 kpc and B1957-I-20, which has V± ~ 173 
km s“^ at a distance of 1.2 kpc (Aldcroft, Romani Sz Cordes 1992). For the most part, these 
objects are consistent with our determination of the 3D rms velocity based on the locations of 22 
MSPs and the absence of MSPs in substantial portions of the volumes searched in high-latitude 
surveys. However, the estimated transverse speed for B1257-1-12 is inconsistent with the overall 
distributions in z and velocity that we have derived, even though it was included in the fitting. 
One possibility is that its distance is overestimated, perhaps by as much as a factor of two, an 
amount sufficient to bring it into consistency with the statistical distribution. It is also possible 
that there are several evolutionary paths for producing MSPs (cf. §p!^), most of which produce 
low-velocity MSPs with others creating rarer, faster MSPs. 

Further study of larger samples of MSP proper motions will result from a combination of 
new surveys, which will discover large numbers of MSPs (cf. §|T^), and use of timing and VLBI 
techniques. Use of the VLBA in conjunction with the Arecibo telescope and the Green Bank 
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Telescope should allow measurement of proper motions for dim and slow MSPs out to a few kpc. 


9. BIRTH RATES OF DISK & HALO MSPS 


For the disk-only model of §^, we have found the column density of MSPs with Pi = 1.56 
ms, Lpi = 1.1 mJy kpc“^ for a plane-parallel model in z to be Nii « kpc“^ (Table ^). The 

implied number of MSPs in a disk of radius with P > Pi and Lp > Lpi is 


NMSpi> -P) > Lp) ~ ^ 10 


Rd 

10 kpc 


P 


1.56 ms 


- 1 ± 0.33 


Lr, 


1.1 mJy kpc 


-2 


- 1 ± 0.2 


(34) 


where the upper and lower values denote the 68% interval. Extrapolation on a per mass basis 
from the local disk surface density to a total disk mass, M^isk: implies 


NMSPi> > Lp) ~ X 10^ 


Lidisk 


4 X W^MqJ V1-56 ms 


P 


- 1 ± 0.33 


Lr. 


1.1 mJy kpc 


-2 


- 1 ± 0.2 


(35) 

These estimates do not include any correction for pulse beaming, whose influence is highly 
uncertain for MSPs. Estimates for this correction range from 1 to 3 (e.g. Bailes & Lorimer 1995). 
The totals are sensitive to the cutoff at small periods and at small luminosities. 


The corresponding birthrate for MSPs, if constant over a galactic age 10^^ yr, for the uniform 


disk 


Nm3p{> U) = i.eiSI X 10-“ y-' (^) . 


and for the extrapolated surface density is 

Nmsp{> Pi) = 3-0;j;2 x 10"® yr"^ 


L^disk 


4 X IOIOM 0 J ' 


(36) 


(37) 


From our constraints on diffuse populations of MSPs, we conclude that, in the vicinity of the 
Sun, the MSP birth rate per unit volume is 100 times less than that from the disk. 


9.1. Comparison with Other MSP Population Studies 

Our estimates may be compared with those derived by Bailes &: Lorimer (1995), Lorimer 
(1995) and Lorimer et al. (1995), who used the Vivekenand & Narayan scale-factor method to 
determine the number of MSPs in the Galaxy and the associated luminosity function. In their 
analyses, specific spatial distributions for the MSPs were adopted to derive the scale factors. BL 
assumed two different scale heights (0.3 and 0.6 kpc) along with a fixed radial distribution to 
estimate lO^’'^ and 10^'® MSPs, respectively, for Lp > 2.5 mJy kpc^ and if all MSPs beam toward 
us. 
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Lorimer et al. (1995) use the radial distribution of Lorimer et al. (1993) (a Gaussian with 
radial scale of 4.8 kpc) and a Maxwellian velocity distribution with rms velocity = \/3 x 100 km 
s“^ to estimate (1.3 ± 0.2) x 10^ MSPs in the Galaxy that are beamed toward us with Lp > 10 
mJy kpc^. 

Lorimer (1995) deduced lower bounds on the scale height and mean 3D space velocity for 
MSPs of 0.5 kpc and 80 km s“^, respectively. These bounds are consistent with our determinations. 

The numbers of pulsars derived by BL and Lorimer et al. (1995) are greater than the estimate 
in Eq. by a factor ~ 2 — 3 for a luminosity cutoff of 2.5 mJy kpc^. Since most of the MSPs 
known are near the Sun (within 2 kpc) an extrapolation to the whole Galaxy is necessary. The 
radial distributions used by BL and Lorimer et al. effectively multiply the uniform disk model 
result by ~ 1.6 and match our own extrapolation (based on scaling up the local disk surface 
density to the given total disk mass in Eq. ^). The extrapolation introduces uncertainty but the 
differences accrue from the following factors. Eirst, the z scale height implied by the Lorimer 
et al. velocity distribution is larger than that derived by us by about a factor of 2. Second, our 
inclusion of scintillation effects yields a search volume that is about 30% larger than otherwise. 
Third, Lorimer et al. include four long period pulsars in their analysis with P > 295 ms that we 
exclude from the MSP sample. Together these differences in the assumed spatial distributions and 
MSP samples explain the size of the differences in the estimated total number of MSPs in the 
Galaxy. 

BL synthesized a luminosity function for MSPs after correcting the observed numbers of 
pulsars for the volume scale factors. Their luminosity function is consistent with a power-law slope 
of —2 (according to our definition of /l^) but with a roll-off below 10 mJy kpc^. Our method is 
able to constrain the lower cutoff on the luminosity function because we evaluate our results at 
values for Lp other than those of actually detected pulsars. 

Similarly, BL suggest that the period distribution decreases in going from 1 to 10 ms and 
roughly estimate that there can be no more than lO'^'^ MSPs with periods with P = 1 ms. Our 
results suggest that the number of pulsars between 1 and 1.5 ms is approximately 50% of the 
number with P > 1.5 ms, or about 5000 pulsars. 


10. RELATIONSHIP TO LOW-MASS X-RAY BINARIES 


10.1. Scale Heights of LMXBs and MSPs 

The evolutionary paths that lead to MSPs are poorly understood (for a review see 


Bhattacharya 1995 ). If all MSPs are “spun up” by mass transfer from a companion star during an 
LMXB phase, then the birth rate of LMXBs must exceed that of MSPs. Kulkarni and Narayan 
(1988) estimated that the birthrate of field LMXBs is about 1-10% of the birthrate of field MSPs 
for an assumed LMXB lifetime of 10® years. With a diminished LMXB lifetime (10^ years), the 
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birthrates are brought into agreement. Our improved estimate of the total number of MSPs in 
the Galaxy does not significantly alter the rate mismatch nor its strong dependence on LMXB 
lifetime. However, our work of the last section does point out that the extrapolation from the 
local MSP population to that of the Galaxy is uncertain by a factor of ~ 2 — 4 (and an additional 
factor of 1 — 3 for beaming). As we will argue below, the kinematic similarity of the LMXB and 
MSP population manifested in the observed scale heights is reasonably strong evidence of an 
evolutionary link; given the great uncertainty in LMXB lifetimes, the best evidence for a causal 
connection between the two populations is not found in the relative number but in the similar 
spatial distribution. In addition, a comparison of the scale height distribution of MSPs with that 
of LMXBs can place significant constraints on evolutionary scenarios leading to MSPs. 

The galactic LMXB scale height derived from analysis of a flux-limited sample ( [Naylor fc 
Podsiadlowski 19^ ) is (0.44-1.17) kpc. Alternatively, based on distance estimates to a subset of 
LMXBs, van Paradijs and White (1995) infer a scale height of 0.5 kpc and, furthermore, argue 
that the LMXBs are predominantly located at Galactocentric radii less than 5 kpc. Although 
the two vertical scales are comparable, the interpretations are quite different. Van Paradijs and 
White assume that the LMXBs have Pop II progenitors and that the scale height is set by large 
velocity kicks at birth (of order 400 km s“^) and the local disk acceleration, which they argue 
is 2.5-4 larger in the relevant inner regions of the Galaxy than locally. This analysis ignores the 
finite birth scale height and the lifetime of the objects. Naylor and Podsiadlowski, on the other 
hand, infer that the LMXBs derive from Pop I stars and have a scale height perpendicular to the 
plane that is roughly like that of the observed thin disk (with a small additional kick), which has 
a nearly constant value. 

Our local determination of the MSP scale height is (0.53-0.81) kpc, comparable with the above 
LMXB estimates. If MSPs are descendants of the LMXB phase and if the scale height increases 
with age, then the local MSP population should have a scale height greater than or equal to the 
LMXB value. If the kicks were as large as suggested by Van Paradijs and White, the minimum 
local scale height of the MSPs would be (1.25-2) kpc, clearly inconsistent with our results. In 
addition, our upper limit on the diffuse number density of pulsars suggests that the observed 
MSPs were born in the disk (Pop I). A consistent interpretation of the LMXB and MSP data is 
that both are Pop I and both are derived from a similar evolutionary channel. 


10.2. Origin of MSP Space Velocities 


In future work, we will discuss how the observed scale height of MSPs and the inferred z 
velocities (~ 50 km s“^) place stringent constraints on the evolution of binary systems that lead 
to MSP formation. One of the main problems in understanding the LMXB evolution path is that 
the formation rate (10“® yr~^) in the Galaxy is so small that the pathway is a priori special. A 
proposed scenario is as follows ( Webbink &: Kalogera 1994] ). A binary composed of a massive 
star (Ml = 10 — 20Mq) and a light companion (M 2 < 0.12Mi) with an initial orbital separation 
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less than about 1000 Rq will pass through an epoch of unstable mass transfer and common 
envelope evolution once the massive star begins to swell. The interaction ejects much of the 
envelope, drawing the pair to very small distances of separation. If the helium core of the primary 
is sufficiently massive it is able to continue to burn and collapse even after its outer hydrogen 
envelope has been removed. The resultant neutron star has an orbit far smaller than the size of 
the original giant primary, an essential requirement if an LMXB phase is to take place. Starting 
with the pre-supernova system we have analyzed how a tight binary is affected by a combination of 
(1) asymmetric SN kick, (2) impact of ejected shell and (3) dissipative processes in the eccentric, 
surviving binary. 

Two effects sculpt the properties of the binaries that survive to give LMXBs and/or MSPs. 
The intrinsic kick given a neutron star by the supernova explosion unbinds loosely bound binaries 
while the impact of the supernova shell on the secondary is responsible for destroying and/or 
unbinding tight binaries. The surviving binaries occupy a relatively narrow range in pre-supernova 
orbital separation, primary and secondary masses and have a limited range of center of mass 
velocities. Most of this analysis is independent of the specific evolutionary pathway leading to the 
pre-SN progenitor. (We present the details of this analysis in Chernoff & Cordes 1996b.) 


11. SELECTION EFFECTS AGAINST FAST BINARIES 


Orbital motion causes MSPs in compact binaries to be missed in surveys that assume the 
pulse period to be constant rather than Doppler shifted (e.g. Johnston &: Kulkarni 19M| ). All 
blind surveys for MSPs, including the 8 analyzed in this paper, make this assumption. One 
circumstance in which the results of previous sections may be altered is if the spin period and/or 
luminosity depend in some way on orbital period. For example, a relation between spin and orbital 


period might be expected on general evolutionary grounds for spun up MSPs (Alpar et al. 1982 


Ruderman Sz Shaham 19^) . Some observations suggest a weak positive correlation ([Lundgren, 


Zepka &: Cordes 1995|) implying that the selection against detecting spin periods less than 1.5 ms 
may be stronger than we have estimated. Because the measured correlation is weak we believe 
that any modification to the distribution of spin periods on this account will be modest. In any 
case, fast binaries have been missed in MSP surveys and their ultimate detection can only increase 
our estimated space densities for MSPs. 

We now give a brief account of survey sensitivity to binary orbital period. The observation 
time T ^ 1 min for most of the Arecibo surveys, so that orbital effects are negligible for 
Port ^ 1.^6 (P in ms) for WD companions with M 2 = 0.3Mq. However, surveys 4,7 & 8 

with T ~ 3 min are insensitive to orbital motion only for ^ 8.^5 Weighted by volumes 

searched, the surveys with longer T contribute strongly to an overall selection against MSP 
binaries with short periods. Indeed, J0751-I-1807 with Porb = 6.^3 was discovered in survey ^ 4 
in a single harmonic, the higher harmonics having been attenuated by orbital motion ( Lundgren] 
Zepka &: Cordes 1995|). That yet-faster binaries with fairly massive WD companions exist is 
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certain because objects like 30751+1807 experience orbital decay due to gravitational radiation 
on less than a Hubble time. Indeed, if MSP-WD binaries achieve Port ^ 8^ solely due to such 
inspiral, then it may be shown that the orbital period distribution dN/dPorb ^ Porb^^^- Overall, 
our conclusions are unaffected for MSPs in binaries with orbital periods ^6 hr and with companion 
masses <0.3Mq. 

Proper consideration of orbital effects — and estimation of the MSP orbital period distribution 
— requires an analysis of search volumes as a function of orbital period and companion mass as 
well as spin period and luminosity. Such a study is beyond the scope of the present paper. 

Future surveys made with greater sensitivities than heretofore using the upgraded Arecibo 
Observatory and the new GET will be able to probe to greater distances while also circumventing 
orbital suppression of Fourier harmonics. Furthermore, algorithms that correct for orbital motion 
are becoming much more feasible with the prospect of computers with teraffops capability. 


12. OPTIMAL SEARCHES FOR MILLISECOND PULSARS 


The population distributions we have derived may be used to optimize new searches for 
MSPs. Search sensitivities (Appendix 0) depend on sky background, dispersion and scattering 
as well as on the flux densities and periods of the MSPs. Consequently, the optimal search is 
frequency and telescope dependent. Here we illustrate the contributions from different effects by 
showing 5Vs, the search volume (the volume searched in a beam area, averaged over P and Lp, 
Eq. 0) and the detection volume, 5Vd, (the volume searched weighted by the dimensionless density 


of MSPs, Eq. 18). To calculate each we use the the best-fit distributions for P and Lp for the 
exponential disk model of and Table By definition SVd < (514: surveys that search much 
more deeply than the scale height of the MSP population yield 6Vd <C (514. For concreteness, 
we consider a survey conducted by telescopes like the under-construction Green Bank Telescope 
(GBT) and a hypothetical analog in the Southern hemisphere, in order that we may consider a 
full sky survey. We assume receiver and survey parameters such that the minimum detectable flux 
density is about 2 mJy when looking at high galactic latitudes and long periods. 


Figure 11 shows (514 and (514 per square degree for a search at 430 MHz. Search volumes 
(top portion of figure) increase more or less monotonically with galactic latitude but level off for 
\b\ > 30°. The volume is smallest toward the galactic center where the sky background is high 
and dispersion and scattering effects are large. By contrast, the detection volume (bottom part of 
figure) is maximum for |fj| ~ 20° and |/| ^ 50° and corresponds to directions that allow the largest 
dj)^ ^ x^ie latitude constraint ensures that the search depth does not exceed the 
MSP scale height. The longitude restriction follows from the variation of the search volume in the 
plane. The detailed shapes of the contours are dependent on the survey frequency and duration 
(per direction) but suggest that future surveys which concentrate on low latitudes will maximize 
the number of new discoveries. However, deep high latitude surveys will better constrain the falloff 
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toward larger z of the disk population as well as place tighter constraints on or make detections of 
any bona fide diffuse or halo-population pulsars. For the hypothetical survey depicted in Figure 
pd| , a total detection volume ~ 13.3 kpc^ is sampled, corresponding to discovery of ~ 585]'^2io 
disk MSPs. This number is for a uniform disk component; any galactocentric radial dependence, 
likely to increase the number of MSPs toward the inner Galaxy, will only increase the number of 
detected MSPs. 


13. DISCOVERING FAST PULSARS 


Our fitting indicates that available survey data already place useful constraints on the 
minimum spin period in the MSP population (cf. §|^ and Table ^). The reason such constraints 
may be placed is found in Figure ^ which shows that, for periods less than 1 ms, a nonzero 
(though small) volume has been searched. Here we estimate how much additional volume must be 
searched in order to expect to find pulsars with Pi < P < Pfast where Pfast is the maximum period 
of interest. 


Using Eq. we derive an upper bound on the volume that must be searched (evaluated at 
the minimum MSP period. Pi), in order that we find pulsars with periods faster than Pfast- Since 
surveys at low periods do not see to large Dmax, we assume that they do not see as far as the z 
scale height ~ 0.5 kpc. Performing the integrals we may write 


{Np) >nd{l 


Pi/Pi: 


’ast ) 




(38) 


Defining the term in square brackets as Vs (Pi) and requiring (Np) ~ 1 to obtain a likely detection, 
we hnd an upper bound on the required search volume to be 


Vs{Pi) 


1 

nd(l - Pl/Pfast) ' 


(39) 


Evaluating Eq. for Pfast = 1-5 ms and Ud ~ 44 kpc“^, we find that, as a function of the 
minimum period, V 5 (Pi) ranges from ~ 1/25 kpc^ for Pf = 0.65 ms (our 99% lower bound on Pf) 
to ~ 1/3 kpc^ for Pf = 1.4 ms. Comparison with Figure shows that, for the luminosity assumed 
for that figure {Lp = 16 mJy kpc^), the Parkes survey (#7) yields search volumes at these Pi that 
are comparable to those needed to yield a detection of a pulsar faster than 1.5 ms. However, the 
assumed Lp for the figure is larger than average and the Parkes survey observes to depths that, 
for some directions, exceed the z scale height. Consequently, it is not surprising that a pulsar 
faster than PSR B1937-I-21 (P = 1.56 ms) has not been found. Nonetheless, future surveys should 
be able to probe this region of period space and either find fast pulsars or determine better the 
period cutoff to the MSP population. Our results indicate that deeper surveys at low galactic 
latitudes (e.g. |5| < 10°) will yield the search volume needed to accomplish these goals. 
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Fig. 11.— (Top:) Aitoff projection in galactic coordinates showing the search volume (kpc^ deg“^) 
for a hypothetical full-sky survey at 430 MHz. The calculated volumes are averages over the period 
and luminosity distributions in our best-fit, exponential disk model. The minimum volume is toward 
the galactic center. The thinnest contour (toward the inner Galaxy) corresponds to the least volume 
( 10 - 3.6 (;| 0 g- 2 ^ while the thickest line (at high latitudes) corresponds to the greatest volume 

(~ 10“^'® kpc^ deg“^). Most of the variation is from a deep minimum toward the Galactic center 
to a shallower variation beginning at | 6 | ~ 15°. The greatest volumes searched are those toward the 
highest latitudes. Structure is seen in the plotted contours (e.g. the North Polar Spur) because the 
search sensitivity and, hence, depth are dependent on the sky background and on dispersion and 
scattering. (Bottom:) Similar projection for the survey detection volume. The minimum contour 
(toward the Galactic center) is 10“^'^ kpc^ deg“^ while the maximum (thickest line) is 10“^'® kpc^ 
deg“^ at latitudes | 6 | ~ 5° and longitudes |£| > 50°. Note that there is no Galactocentric radial 
dependence of our assumed MSP number density, so all the structure at 6 = 0° (and other lines of 
constant latitude) is due to the depth of the survey. 
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14. DISCUSSION 

Through a likelihood analysis, we have constrained the period and pseudo-luminosity 
distributions to be steep power laws with slopes ~ —2. The distributions imply that the 
population of MSPs increases rapidly to smaller periods and smaller pseudo-luminosities. We 
infer a minimum period Pmin > 0.65 ms at 99% confidence and a minimum luminosity cutoff 
Lpi = mJy kpc^. The column density of MSPs in the local vicinity of the solar system is 

Nd ~ 50^20 kpc~^- The limits on a diffuse halo-like component are ^1% of the midplane density. 
All these results are essentially identical for each of the models we have analyzed. Estimates of the 
total number of MSPs in the Galaxy are uncertain. Extrapolating on a per mass basis from the local 
disk surface density to a total disk mass, M^isk, we find Nmsp ~ ^ 10^ [Mdisk/4: x 

for cutoff period 1.56 ms and cutoff luminosity 1.1 mJy kpc“^ (without correction for beaming). 

Our analysis assumes specific forms for the period and luminosity distributions, namely 
power-law functions, that undoubtedly influence the specific values for numbers of pulsars in a 
given period range and also on the minimum period cutoff. We have not tested other mathematical 
forms for these distributions, so the true cutoff for the period distribution may be different than we 
have derived. Nonetheless, because the period distribution montonically increases with decreasing 
period, our quoted minimum period is larger than it would be for a function that plateaus or 
decreases with decreasing period below 1.56 ms. We consider the most important implication of 
our derived minimum period to be that MSPs faster than those already found may indeed be 
present in the Galaxy: the surveys done heretofore cannot rule out their existence. In addition, 
modeling of the spinup process using full general relativity (Gook et al. 1994a,b) implies that 
gravitational instabilities do not prohibit the formation of very fast MSPs. Of course the ultimate 
existence proof for MSPs with P < 1.56 ms lies in future surveys that can explore large volumes 
of the Galaxy at these small periods. Such surveys will be feasible with new spectrometers that 
can sample more frequency channels at faster rates and with post processing that can contend 
with motion of fast pulsars in binaries. 

Another implication of our results on the period distribution is that, if MSPs exist due to 
accretion-driven spinup of neutron stars, then accretion must ensue for sufficiently long times 
that periods shorter than 1.56 ms can be achieved. From the work of Gook et al. (1994a,b), 
such accretion appears possible without requiring typical ages for LMXBs that are so long as to 
resurrect the discrepancy between birth rates for MSPs and LMXBs. 

The observed scale height of MSPs implies that they are a low-velocity population among 
neutron stars, having an rms speed that is about a factor of 5 smaller than that of young pulsars 
with much stronger magnetic fields. A part of the total inferred dispersion we attribute to a kick 
unique to the evolution of MSP systems (~ 40 km s“^) and the rest to the effect of diffusive 
processes that increase the dispersion of old objects. A number of kinematic signatures that 
should be evident in larger MSP samples (transverse motions, asymmetric drift, shape of velocity 
ellipsoid) are described. 
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The disparity in velocity between the low-field MSPs and high-held pulsars might be taken as 
evidence that the two empirical classes of neutron stars are born through substantially different 
processes. If MSPs are produced largely through accretion-induced collapse of a white dwarf and 
if that process yields only a small kick to the resultant NS compared to Type II supernova then 
the observed dispersions of MSPs and high-held pulsars may hnd a natural explanation. In any 
case, the similarity in scale height of MSPs and LMXBs shows that the formation of the NS in 
both objects is accomplished without substantial center-of-mass impulses and supports the notion 
of an evolutionary connection. On the other hand, if binary survival after the type II supernova 
is the most signihcant bottleneck in the production of LMXBs and their MSP descendants, it 
is possible that the processes that dictate survival of the binary system are also responsible for 
allowing only a limited range of center-of-mass velocities. Correlations between spin and orbital 
periods and space velocity, such as those suggested by Bailes et al. (1994), depend critically on 
the details of mass transfer and on the number of evolutionary paths that lead to MSP formation. 
Elsewhere, we will present our detailed analysis of the effects that sculpt binary survival. 
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APPENDICES 


A. SEARCH SENSITIVITIES 


The pulsar searches we consider involve the removal of dispersion delays between the outputs 
of a multichannel receive using trial values for the dispersion measure. The resultant time series 
is then Fourier analyzed. Suppose an Nppx-length Fast Fourier Transform (FFT) is calculated 
from the time series for each trial dispersion measure. With a sample time At and pulse period P, 
harmonics appear in frequency bins 


ke = 


i^tNppp 

P ’ 


£ = 0 , 1 ,... 


(Al) 


including a “DC” term (£ = 0) and the fundamental (£ = 1). Let the intrinsic pulse shape be 
s{t), 0 < t < P so that a (short) discrete Fourier transform (DFT) of this shape over a single 
pulse period is s{i), i = 1,... ,M. This function would determine the envelope of harmonic 
amplitudes in the long FFT were it not for additional contributions that derive from the 
post-detection averaging time (or “time constant”), from dispersion smearing across individual 
frequency channels, and from pulse broadening due to interstellar scattering (for distant sources). 
In many surveys, post-detection smoothing is simply an RC filter whose time domain response is 
a one-sided exponential function. Interstellar scattering produces nearly the same kind of time 
response, while the dispersion time function is dictated by the shapes of receiver filters, usually 
approximately Gaussian in form. Letting the M-point DFTs of the time constant, dispersion, and 
scattering functions be stc, Sd and Sg, respectively, we may write the effective envelope function of 
harmonics as 

^effif-) = s{£)stc{i)Sdii)Ssii)- (A2) 

It is useful to define the ratio of the £-th harmonic to the DC value as 


Ri 


Seff{£) 

^e//(0) 


(A3) 


Survey FFTs are analyzed by constructing partial sums of harmonics (of the FFT magnitude 
or squared magnitude) for different trial periods. These sums are typically of Nh = 1,2,3,4 ,8 and 
16 harmonics, though there are variations on this. Suppose that a threshold rjp is chosen that 
represents the number of standard deviations in the FFT’s magnitude. This is typically 77 ^ ~ 6 
to 9 in order to minimize false-alarms when testing large numbers of spectral values (typically 
multiples of 10 ®) in a survey. 
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The minimum detectable flux density for a sum of harmonics 1,... , Nh is 

Cl _ __ 

min,Nh NpolAl^ AtNpFT 

\£=1 / 

(where Tgyg is the system temperature [K]; G is the telescope gain [K Jy“^]; Npoi = 2 is the number 
of independent polarization channels included; Av is the total bandwidth; and At is the sample 
interval). For searches that analyze |FFTp rather than |FFT|, Ri —> i?|. The actual minimum flux 
density depends on the number of harmonics that contribute significantly which, in turn, depends 
on the duty cycle of the pulse. Because extrinsic effects (viz. dispersion and scattering), broaden 
the pulse, the optimal and corresponding Smin are strongly dependent on the observation 
frequency, distance, direction and pulse period. The direction dependence is manifested in the 
dispersion measure to which a pulsar of given period may be detected. Consequently Smin is the 
minimum over all considered in the analysis and may be written with dependences 


V7^ 

Nh 

TRi 


(A4) 


Smin — Smin {£, b, P, DM, v, Av, Nch 1 Tgys , G, • • • ). (A5) 

In practice, surveys usually test only a subset of all possible harmonic sums. We take this into 
account when computing Smin for each survey. 

Note that our expression for the minimum flux density differs from that often quoted in the 

with a factor yJw/{P — w), where w is the pulse width. The divergence of this factor as w ^ P 
is equivalent to assuming Ri = 0, which overestimates the true Smin because, even when pulse 
smearing exceeds a pulse period, the variable flux remaining at the fundamental frequency can 
still be detectable for a luminous pulsar. Our expression takes this possibility into account, which 
corresponds to 0 < Pi <C 1. 

It is important to calculate accurately the minimum detectable flux density because it 
determines the galactic volume searched. This volume is small but not zero for very short periods 
<1.5 ms. 


literature (e.g. Camilo, Nice &: Taylor 1996|), which replaces the factor in large brackets in Eq. A4 


B. INTERSTELLAR SCINTILLATIONS 

Interstellar scintillations are intensity variations in both time and frequency caused by 
multipath propagation through ionized gas. At 400 MHz, both diffractive (DISS) and refractive 
(RISS) interstellar scintillations contribute to the flux variations of pulsars. Here we restrict the 
discussion to DISS, which will dominate RISS at 400 MHz and is especially important because its 
probability density is skewed whereas RISS is symmetric. 
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Log Luminosity L^ (mJy kpc^) 

Fig. Bl.— Luminosity functions with and without the effects of scintillations included. (Heavy 
Solid Line:) Intrinsic luminosity function having a power-law slope of —2. (Light Solid Line:) The 
scintillated luminosity function when interstellar scintillations are saturated and unquenched by 
time-bandwidth averaging. (Dotted Lines:) Scintillated luminosity functions with various degrees 
of averaging, indicated by riiss (cf. Eq. |B^) . 
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DISS causes the pulsar flux density S to vary as S' = gS, where g is the DISS gain that has a 
one-sided exponential distribution when DISS is saturated and not quenched by time-bandwidth 
averaging (Rickett 1990; Cordes &: Lazio 1991). Except for the very nearest pulsars {D < 100 pc), 
DISS is saturated at 400 MHz. The characteristic time and frequency scales of DISS diminish with 
increasing distance. Use of finite bandwidth B and data-span length T will average over distinct 
scintillation maxima, increasing the number of degrees of freedom from 2 (for unquenched DISS) 
to 2 niss, where 

n,,, ~ (l + 0.2^) (^1 + 0.2^) , (Bl) 

is the characteristic bandwidth of DISS, and is the characteristic time scale (Cordes 
1986). The characteristic bandwidth and time scale have been measured for many pulsars and 
were used as input to the Taylor &: Cordes (1993) model for pulsar distances. For our purposes, 
we use the TC model’s estimation of the scattering measure along with the distance and frequency 
to estimate the scintillation parameters. 

The pdf of g is 

fgig^n^ss) = (B2) 

SE [riiss) 

with U{g) the Heaviside function and T is the gamma function. As riiss —> oo (i.e. pulsars at large 
distances or observed at low frequencies), fg tends toward a delta function, 5{g — 1). 

We include scintillations in our analysis by defining a scintillated pseudo luminosity, L'p = gLp. 
For a luminosity function fLp{Lp), the corresponding scintillated luminosity function is 

fVpiL'p) = j dgg-^fg{g,niss)fLp{L'p/g). (B3) 

The distinctive effect of DISS is that, if the intrinsic luminosity function has cutoffs at low and 
high luminosities, the scintillated luminosity function will not. In fact, the scintillated luminosity 
function will extend to zero luminosity because the most probable scintillation gain (for riiss = 1 ) 
is zero. Luminosities larger than the upper cutoff will be seen owing to the long exponential 
tail of fg (again for riiss = !)• Figure shows examples of scintillated luminosity functions for 
several values of riiss- As riiss oo, the scintillated luminosity function tends toward the original, 
unscintillated luminosity function. 

In surveys where DISS is saturated and unquenched {njss = 1) and single trials are made on 
each sky position, the volume surveyed is effectively increased by a factor ( 5 ^/^) = r(|) ~ 1.33. 
Multiple trials can increase or decrease this volume factor, depending on how the results of the 
various trials are combined. 
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Table 1; MSP Survey Parameters 


Survey 

Site 

V 

VLs 

Nmsp 

Ssys 

*5mino 

Ref 



(GHz) 

(deg2) 


(Jy) 

(mJy) 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

1 

A 

0.43 

680 

3 

3 

0.5 

1 

2 

A 

0.43 

235 

4 

3 

1.0 

2 

3 

A 

0.43 

250 

0 

3 

0.4 

3 

4 

A 

0.43 

7 

1 

3 

0.2 

4 

5 

A 

0.43 

682 

2 

3 

0.7 

5 

6 

A 

0.43 

150 

1 

3 

0.4 

6 

7 

P 

0.44 

20,600 

10 

90 

3.0 

7 

8 

J 

0.41 

1,650 

1 

70 

3.1 

8 


Sites: A = Arecibo, J = Jodrell Bank, P = Parkes. 

References: (1) Camilo et al. 1996; (2) Nice et al. 1995; (3) Thorsett et al. 1993; (4) Lundgren et 
al. 1995; (5) Foster et al. 1995; (6) Wolszczan 1990; (7) Manchester et al. 1996; (8) Nicastro et al. 
1995. 
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Table 2; Millisecond Pulsars Used 


MSP Name 

(1) 

1 

(deg) 

(2) 

b 

(deg) 

(3) 

P 

(ms) 

(4) 

log AP 
(ms) 

(5) 

S 

(mJy) 

(6) 

AP 

(mJy) 

(7) 

Dl 

(kpc) 

(8) 

Du 

(kpc) 

(9) 

Ref 

(10) 

J0034-0534 

111.5 

-68.1 

1.88 

-10.7 

16 

5 

0.74 

1.23 

7 

J0437-4715 

253.4 

-42.0 

5.76 

-11.4 

600 

180 

0.105 

0.175 

7 

J0613-0200 

210.4 

-9.3 

2.19 

-11.4 

21 

6 

1.64 

2.74 

7 

J0711-6830 

279.5 

-23.4 

5.49 

-4.1 

7 

2 

0.77 

1.29 

7 

J0751+1807 

202.7 

21.1 

3.48 

-11.0 

10 

3 

1.51 

2.53 

4 

J1012+5307 

160.3 

50.9 

5.26 

-10.7 

30 

9 

0.39 

0.65 

8 

J1045-4509 

280.9 

12.3 

7.47 

-10.7 

20 

6 

2.43 

4.05 

7 

B1257+12 

311.3 

75.4 

6.22 

-12.7 

20 

6 

0.47 

0.77 

6 

J1455-3330 

330.7 

22.6 

7.99 

-10.2 

13 

4 

0.56 

0.93 

7 

J1640+2224 

41.1 

38.3 

3.15 

-12.3 

12 

4 

0.88 

1.48 

5 

J1643-1224 

5.7 

21.2 

4.62 

-10.5 

75 

23 

4.84 

oo 

7 

J1713+0747 

28.8 

25.2 

4.57 

-12.1 

36 

10 

0.8 

1.6 

5 

J1730-2304 

3.1 

6.0 

8.12 

-10.5 

43 

13 

0.38 

0.64 

7 

B1855+09 

42.3 

3.1 

5.36 

-12.5 

31 

9 

0.70 

1.30 

2 

B1937+21 

57.5 

-0.3 

1.56 

-12.7 

240 

72 

3.60 

15.7 

2 

B1957+20 

59.2 

-4.7 

1.61 

-12.52 

20 

6 

1.15 

1.91 

2 

J2019+2425 

64.7 

-6.6 

3.93 

-12.7 

15 

5.0 

0.68 

1.14 

2 

J2124-3358 

10.9 

-45.4 

4.93 

-10.2 

20 

6 

0.18 

0.30 

7 

J2145-0750 

47.8 

-42.1 

16.05 

-9.7 

50 

15 

0.38 

0.62 

7 

J2317+1439 

91.4 

-42.4 

3.45 

-12.7 

14 

5 

1.4 

2.3 

1 

J2322+2057 

96.5 

-37.3 

4.81 

-12.6 

4 

2 

0.5 

1.1 

1 

J2229+2643 

87.7 

-26.3 

2.98 

-12.1 

18 

5 

1.0 

2.0 

1 


References: (1) Camilo et al. 1996; (2) Nice et al. 1995; (3) Foster et al. 1995; (4) Lundgren et al. 
1995; (5) Thorsett et al. 1993; (6) Wolszczan 1990; (7) Manchester et al. 1996; (8) Nicastro et al. 
1995. 
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Table 3; 


Best-fit Disk Models 


Parameter 

Gaussian 

in z 

Exponential 

in z 

Gaussian 

in V 

units 

Otp 

2.0± 0.33 

2.0±0.33 

2.0±0.33 

— 


2.0± 0.2 

2.0±0.2 

2.1±0.2 

— 

CTz 

0.65^0:16 

o.5ol°:ll 

— 

kpc 

^z,birth^ 

— 

— 

0.1 

kpc 

av 

— 

— 

52111 

km s“^ 

nd 

29^11 


IT q+28 

00-18 

kpc“^ 

Nd 

iro+29 

0Z_19 


49127 

kpc“2 

Pi 

> 1.0 (95%) 

> 0.65 (99%) 

> 1.0 (95%) 

> 0.65 (99%) 

> 1.0 (95%) 

> 0.70 (99%) 

ms 

Lp^ 

1 1+°-^ 

1 1+°-^ 

-‘-•-‘--0.5 

1 1+°-^ 

-‘-•-‘--0.5 

mJy kpc2 


Confidence intervals, except where noted, are two-sided 68% intervals. ^ fixed parameter. 
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Table 4; 


Disk + Diffuse Models 


Model 

Diffuse only 


Disk + Diffuse 


nh 

(kpc“^) 

£(0,%) 

£(nrf,0) 

nd 

(kpc“^) 

nh 

(kpc“^) 

nh{Ro) 

(kpc“^) 

uniform 

1.5 

10 - 21.3 

38 

< 0.42 (90%) 


density 




< 0.84 (99%) 


rh = l kpc 

4080 

10 - 16.4 

38 

< 1520 (90%) 

< 0.83 

Sh = 7/2 




< 3040 (99%) 

< 1.66 

r/j = 5 kpc 

10.8 

10-15.1 

38 

< 3.3 (90%) 

< 0.31 

Sh = 2 




< 6.7 (99%) 

< 0.62 
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Table 5; Velocity Moments for Nearby Pulsars 










ay 

0-2 

Vr 

^VR 

Vt 

6vt 

Vz 

SVz, 

(km s“^) 

(kpc) 







Uniform Surface Density (independent of R): 



20 

0.05 

0.02 

1.28 

0.00 

0.86 

-0.00 

0.70 


0.15 

0.01 

1.28 

-0.00 

0.86 

-0.00 

0.79 

40 

0.05 

0.00 

1.26 

-0.11 

0.80 

-0.00 

0.66 


0.15 

0.01 

1.26 

-0.11 

0.80 

-0.01 

0.69 

60 

0.05 

0.00 

1.18 

-0.22 

0.77 

0.00 

0.64 


0.15 

0.00 

1.18 

-0.22 

0.77 

-0.00 

0.65 

80 

0.05 

0.01 

1.09 

-0.29 

0.75 

0.00 

0.61 


0.15 

0.00 

1.09 

-0.30 

0.75 

0.00 

0.61 

100 

0.05 

0.01 

1.01 

-0.35 

0.74 

0.00 

0.58 


0.15 

0.00 

1.00 

-0.35 

0.74 

0.00 

0.58 

Exponential Surface Density (in 

R)-. 




20 

0.05 

0.03 

1.24 

-0.17 

0.83 

-0.01 

0.69 


0.15 

0.01 

1.23 

-0.17 

0.83 

-0.01 

0.79 

40 

0.05 

-0.01 

1.17 

-0.31 

0.78 

-0.00 

0.63 


0.15 

0.01 

1.17 

-0.31 

0.78 

-0.00 

0.66 

60 

0.05 

-0.01 

1.09 

-0.41 

0.73 

-0.00 

0.58 


0.15 

0.00 

1.09 

-0.41 

0.74 

0.00 

0.60 

80 

0.05 

0.00 

1.03 

-0.47 

0.69 

-0.00 

0.55 


0.15 

0.01 

1.02 

-0.47 

0.70 

0.00 

0.56 

100 

0.05 

0.01 

0.97 

-0.50 

0.66 

-0.00 

0.51 


0.15 

0.00 

0.96 

-0.50 

0.66 

0.00 

0.52 


Velocity moments for particles of Galactocentric radius 7.5 < R < 9.5 kpc and jz] < 3 kpc. Here 
vr means < vr > /ay, Svr means \/< {vr— < vr >)2 >fav, and so on. The top section refers 
to a disk with constant birth density in R; the bottom section to a disk with birth density varying 
with exponential scale length 3.5 kpc in R. All moments given in units of ay- Numerical accuracy 
of ±0.02 for all entries. All mixed second order moments are zero to ±0.03. 



